## summary of the output from the 1d model with alpha = 1
## (balanced data)
##
##

library(coda)
load("pairData-AttnCheckPassed.Rda")


load("out-1d-fixedalpha-nofraud-balanced-newIDs-AttnCheckPassed.Rda")
out <- out.1d.fixedalpha

thetas <- out[, grep("theta", colnames(out))]
## alpha <- out[, grep("alpha", colnames(out))] not necessary (all == 1)

M <- nrow(out) ## number of MCMC draws

## calculate respondent-statement-specific in-sample prediction error rates
pairData.aug <- pairData
pairData.aug$post.prob.correct <- NA
for (i in 1:nrow(pairData.aug)){
    theta.vec.s1 <- thetas[,paste("theta",
                                  pairData.aug$statementID1[i], sep=".")]
    theta.vec.s2 <- thetas[,paste("theta",
                                  pairData.aug$statementID2[i], sep=".")]

    if (pairData.aug$choice[i] == pairData.aug$statementID1[i]){
        pairData.aug$post.prob.correct[i] = mean(pnorm(theta.vec.s1 -
                                                       theta.vec.s2))
    }
    if (pairData.aug$choice[i] == pairData.aug$statementID2[i]){
        pairData.aug$post.prob.correct[i] = mean(pnorm(theta.vec.s2 -
                                                       theta.vec.s1))
    }

    if (i %% 1000 == 0){
        cat("i =", i, "  of", nrow(pairData.aug), "\n")
    }
}




## aggregate to respondent
resp.id <- sort(unique(pairData.aug$respondent))
n.resp <- length(resp.id)
resp.post.prob.correct <- rep(NA, n.resp)

for (i in 1:n.resp){
    data.sub <- pairData.aug[pairData.aug$respondent == resp.id[i],]
    resp.post.prob.correct[i] <- mean(data.sub$post.prob.correct)
}

hist(resp.post.prob.correct, nclass=100, col=rgb(0, 0.4, 0.7))



## aggregate to statement
statement.id <- sort(unique(c(pairData.aug$statementID1,
                              pairData.aug$statementID2)))
n.statements <- length(statement.id)
statement.post.prob.correct <- rep(NA, n.statements)

for (i in 1:n.statements){
    data.sub <- pairData.aug[pairData.aug$statementID1 == statement.id[i] |
                            pairData.aug$statementID2 == statement.id[i],]
    statement.post.prob.correct[i] <- mean(data.sub$post.prob.correct)
}



hist(statement.post.prob.correct, nclass=20, col="dodgerblue")





## compare rank order correlation of theta point estimates and objective truth
## to rank order correlation of theta point estimates and valence
##


statements <- read.csv("./selectedStatementsPunctuationFixed-balanced.csv", header=TRUE, stringsAsFactors=FALSE)

statements <- statements[statements$keep==1,]

valence <- statements$valence
valence[valence=="R"] <- 1
valence[valence=="N"] <- 0
valence[valence=="L"] <- -1
valence <- as.numeric(valence)

obj.truth.6pt <- rep(NA, nrow(statements))
obj.truth.6pt[statements$rating == "pants-fire"] <- 0
obj.truth.6pt[statements$rating == "FALSE"] <- 1
obj.truth.6pt[statements$rating == "barely-true"] <- 2
obj.truth.6pt[statements$rating == "half-true"] <- 3
obj.truth.6pt[statements$rating == "mostly-true"] <- 4
obj.truth.6pt[statements$rating == "TRUE"] <- 5


theta.hat <- colMeans(thetas)


print(cor(valence, theta.hat, method="spearman"))
print(cor(obj.truth.6pt, theta.hat, method="spearman"))


